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Diffusion Monte Carlo (DMC) calculations are performed on the monocyclic and bicyclic forms 
of m-benzyne, which are the equilibrium structures at the CCSD(T) and CCSD levels of coupled 
cluster theory. We employed multi-configuration self-consistent field trial wave functions which are 
constructed from a carefully selected 8-electrons-in-8-orbitals complete active space [CAS(8,8)], with 
CSF coefficients that are reoptimized in the presence of a Jastrow factor. The DMC calculations 
show that the monocyclic structure is lower in energy than the bicyclic structure by 1.9(2) kcal/mole, 
in excellent agreement with the best coupled cluster results. 

PACS numbers: 02.70.Ss, 71.15.-m, 31.25.-v 



I. INTRODUCTION 

Quantum chemistry methods used to study 1,3- 
Didehydrobenzene, commonly known as meta-benzyne or 
m-benzyne, do not agree upon its geometrical structure. 
The issue at question is, which of the structures shown 
in Fig. I, the monocyclic (a) or the bicyclic (b), is lower 
in energy. These two structures differ mostly in the dis- 
tance between the radical centers, CI and C3; the single 
occupied orbitals on carbon atoms 1 and 3 form a bond 
in (b) which decreases the biradical nature of (a). 

At the time of its synthesis coupled cluster calcu- 
lations with single and double excitations plus a pertur- 
bative treatment of triple excitations CCSD(T), with a 
small 6-31G(d,p) basis were used in identifying the 
structure of m-benzyne. The computed infrared frequen- 
cies of the monocyclic structure, (a) , were in good agree- 
ment with the experimental one, thus indicating that this 
structure is the stable one. Subsequent calculations using 
larger basis sets also supported this conclusion [1, 0, [j| . 

In an extensive recent study, Smith, Crawford, and 
Cremer [j| performed coupled cluster CCSD (single and 
double excitations), CCSD(T), and CCSDT (single, dou- 
ble and triple excitations) calculations on m-benzync. 
The CCSD and CCSD(T) calculations were done with 
6-31G(d,p) [|J and cc-pVTZ @ basis sets with refer- 
ence determinants from restricted Hartree-Fock (RHF), 
unrestricted Hartree-Fock (UHF) and Bruckner orbitals. 
They found that CCSD favors the bicyclic form (b) as 
the ground state with a C1C3 bond length sa 1.56 A. 
On the other hand, CCSD(T) predicted (a) as the sta- 
ble structure with a C1C3 bond length w 2.10 A. With 
the small 6-31G(d,p) basis set and with RHF orbitals, 
CCSDT also favored the structure (a). One of the con- 
clusions of their study is that the inclusion of triple exci- 
tations (as in CCSD(T) or CCSDT) is essential to obtain 
a correct energy ordering of the structures of m-benzync, 
and, inclusion of d- and /-type polarization functions in 
the basis is needed for a quantitatively correct energy 
difference. 

The structural predictions of density functional the- 
ory (DFT) are inconclusive, depending strongly on the 
choice of the approximate exchange-correlation func- 
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FIG. 1: Proposed structures of m-benzyne: (a) the mono- 
cyclic structure and (b) the bicyclic structure. 



tional. The pure gradient corrected functionals, Becke- 
Lee- Yang-Parr (BLYP) @, H and Becke-Perdew-Wang- 
91 (BPW91) @, [|[ predict that the monocyclic form, 
(a), is stable, whereas, the more complicated 3-paramcter 
B3LYP functional [HI, which includes an admixture of 
exact exchange, favors the bicyclic structure, (b) [H, 0, 
[Tl| . Since Hartree-Fock theory strongly favors the bi- 
cyclic structure, (b) 0], it is not surprising that the hy- 
brid functional would also. Further, it was found that 
the infrared frequencies computed with the B3LYP func- 
tional of the two structures agree about equally well [ll[ 
with the experimental ones, which led the author to ques- 
tion the previous interpretation [l[ of the experimental 
results. 

Both coupled cluster and DFT energies are based on 
single-reference wave functions and thus their predictions 
are questionable if m-benzyne has considerable multi- 
reference character. Multi-configuration self-consistent 
field (MCSCF) calculations, which do a good job of de- 
scribing static correlation but a poor job on dynamic 
correlation, indicate that (a) is more stable by a large 
amount, 12.9 kcal/mole This is to be contrasted with 
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TABLE I: Optimized geometrical parameters of structures 
(a) and (b) with a cc-pVTZ basis set from Ref. The first 
geometry was optimized using CCSD(T) and the second one 
with CCSD. Bond distances are in angstroms and angles are 
in degrees. 
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the RHF theory which predicts that structure (bHs more 
stable than structure (a) by about 18 kcal/mole Q- This 
shows the importance of a multi-determinant description 
of m-benzyne; and, in order to get reliable energy dif- 
ferences it is essential to have an accurate description of 
both static and dynamic correlations. 

In this paper we use variational Monte Carlo (VMC) 
and diffusion Monte Carlo (DMC) to study this issue. 
These methods have the advantage that both static and 
dynamic correlation can be described accurately, pro- 
vided that well optimized, multi-determinant trial wave 
functions are employed. Further, the results have only a 
weak dependence on the size of the basis. Currently, 
it is still expensive to optimize the geometry of large 
molecules using QMC methods mainly because of their 
stochastic nature. Thus, we limited our study to the 
two geometries of (a) and (b) which are respectively the 
equilibrium geometries for CCSD(T) and CCSD with a 
cc-pVTZ basis @. The geometrical parameters are re- 
produced in Table Q] for convenience. 



II. QMC METHODS 

We employ two QMC methods, variational Monte 
Carlo (VMC) and diffusion Monte Carlo (DMC) which 
have been reviewed before [H, EH, EH ■ In both methods 
the expectation values are obtained by a Monte Carlo in- 
tegration and consequently there is great freedom in the 
choice of the trial wave function. VMC yields expecta- 
tion values for the trial wave function, whereas in DMC 
the trial wave function serves as a starting point for a 
stochastic projection onto the ground state. Although 
DMC is exact for Bosonic ground states, for Fermionic 
systems the projection is performed while imposing the 
boundary condition that the nodes of the many-body 
wave function are the same as those of the trial states. 
Consequently, the DMC energy is also an upper bound 
to the true energy and the fixed-node error depends cru- 



cially on the quality of the trial wave function. In prac- 
tice, we use an accelerated Metropolis algorithm with 
a very small auto-correlation time [15[ and a DMC algo- 
rithm [l6j , that takes into account the singularities in the 
time-evolution operator and allows one to use large time 
steps while still having an acceptably small time-step er- 
ror. All the QMC calculations are performed using the 
program package CHAMP [I?) • 

The trial wave function is of the standard Jastrow- 
Slater form consisting of a Jastrow factor J(R) = e-^ R ) 
times a linear combination of determinants. Specifically, 

* t (R) = /(R)^qC,(R) (1) 

i 

where R = (ri , T2, ■ ■ ■ , r/v) is a 3N dimensional vector, 
is the position vector of the kth electron, Cj(R) is the ith 
configuration state function (CSF), which is a symmetry- 
adapted linear combination of Slater determinants, and 
d are the CSF linear expansion coefficients. The Jas- 
trow part includes explicitly electron-nucleus, electron- 
electron, and electron-electron-nucleus correlation func- 
tions |18|. The statistical and the systematic errors in 
the VMC and DMC energies depend on the quality of the 
trial wave function. The Jastrow and CSF variational pa- 
rameters are optimized using recently developed energy 
minimization methods [l|| Ho, Hl[ . 

To eliminate the core electrons, we employ the energy- 
consistent scalar-relativistic Hartree-Fock pseudopoten- 
tials of Burkatzki, Filippi and Dolg [l^j and the accom- 
panying basis sets. For the DMC calculations, we use the 
standard locality approximation which [23j introduces 
another small systematic error. 



III. CONVERGENCE TESTS 

Before performing the final calculations, we studied the 
convergence of the energy with respect to a) the size of 
the basis, b) the cutoff used for selecting which CSFs are 
included in the expansion, and c) the Trotter time-step 
used to evolve the system. All of our DFT and MCSCF 
calculations are performed using GAMESS [24j ]. 

1. Basis set convergence 

The VMC and DMC energies depend on the quality 
of the basis used to expand the orbitals in the determi- 
nants. One of the advantages of both QMC methods, 
particularly DMC, is that the dependence on the basis 
size is much weaker than for other quantum chemistry 
methods. 

Figure [5] shows the QMC energies obtained, with the 
determinant part of the wave function chosen to be a 
B3LYP single determinant, for various basis sets. DZ 
denotes a double-zeta basis set, TZ* a triple-zeta basis 
set without the /-functions of the carbon atoms, TZ is a 
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FIG. 2: VMC and DMC energies for different basis sets ob- 
tained with a B3LYP trial wave function for the structures 
(a) and (b). The quality of basis sets range from double-zeta 
(DZ) to quadruple-zeta (QZ) for both carbon and hydrogen. 
TZ* is the triple-zeta basis set without the carbon /-functions, 
and QZ* is the quadruple-zeta basis set without the carbon 
(/-functions. The inset shows AE — Et — E a in mHa. 



triple-zeta basis, and QZ* denotes a quadruple-zeta basis 
without the carbon ^-functions. As expected, the VMC 
and DMC energies improve as the basis is improved. The 
DMC energies obtained from DZ and TZ bases differ by 
about 8 milli Hartrce (mHa) and those obtained from 
TZ and QZ* bases by less than 2 mHa. The w 4 mHa 
difference in DMC energy between the TZ and TZ* basis 
sets shows that it helps to include the polarization f- 
functions. 

The quantity of interest is the difference in the energies 
of the two structures and this converges faster than the 
individual total energies. The inset of Fig. [2] shows the 
energy difference AE = Eb — E a is already converged 
within statistical errors even with the double-zeta basis. 
In the rest of our QMC calculations, we employ only 
double- and triple-zeta basis sets. 



2. CSF convergence 

As seen from Figure [5J when a B3LYP single- 
determinant is employed in the wave function, the energy 
of the (b) structure is lower than that of the (a) structure. 
However, since this is a system where static correlation is 
important, it is imperative to employ multi-determinant 
wave functions in QMC. In fact inclusion of additional 
determinants results in reversing the ordering of the en- 
ergies. To this end, we employed truncated complete ac- 
tive space self consistent field (CASSCF) wave functions 
to construct the determinant part of our trial wave func- 
tions, which yield lower VMC and DMC energies than 
the single-determinant wave functions. The active space 
was carefully selected such that it contains the six it elec- 



trons and the two unpaired electrons in structure (a) that 
form the C1C3 bridging bond in structure (b). 

In order to reduce the number of determinants in the 
QMC calculations, we include only CSFs which have a 
coefficient larger in absolute magnitude than a cutoff pa- 
rameter, rj, in the MCSCF wave function. The CSF co- 
efficients are then reoptimized within VMC in the pres- 
ence of the Jastrow; typically, all the CSF coefficients 
(except the dominant one) get smaller in absolute mag- 
nitude upon rcoptimization because of the presence of the 
Jastrow factor. Table HI1 shows the QMC energies for sev- 
eral values of the cutoff rj. With a cutoff of rj = 0.01, the 
DMC total energies are converged to better than 1 mHa, 
and this is the value of r) we used in the rest of our calcu- 
lations. Moreover, there is a good degree of cancellation 
of error between the energies of the two structures, so 
that the energy differences are converged, within statis- 
tical errors, even with larger values of r/. 



3. Time-step convergence 

The imaginary-time evolution operator used in DMC 
is exact only in the limit of zero Trotter time step, St. 
Consequently, energies evaluated at finite St have a time- 
step error which may be positive or negative. Table IIIII 
shows the DMC energies for the two structures and their 
difference for various values of St. The time-step errors 
of the energies of the two structures cancel out within 
statistical error upon taking the energy difference. Con- 
sequently, we employed a time step of St = 0.1 Hartree -1 
in the rest of our calculations. 



IV. RESULTS 

Table IIVI shows the energy difference of the two struc- 
tures AE = E b -E a obtained from the BLYP and B3LYP 
density functional methods and from MCSCF/CAS(8,8) 
theory. The density functional calculations are con- 
verged with respect to basis size; the difference between 
the triple- and quadruple-zeta basis sets being less than 
0.1 kcal/mole. In the CAS(8,8) calculations, we used 
double- and triple-zeta basis sets. In agreement with pre- 
vious papers, B3LYP favors structure (b), BLYP favors 
structure (a), and MCSCF greatly favors structure (a). 

The energy differences AE at the B3LYP and BLYP 
levels of theory are in good agreement with the values 
of —2.9 and 3.0 kcal/mole obtained using an all electron 
6-311++G(3df,3pd) basis set 0]. Also, the CAS(8,8) en- 
ergy difference with the triple-zeta basis is in very good 
agreement with the all-electron value, 12.9 kcal/mole, ob- 
tained using a 6-31+G(2df,2p) basis value at the same 
level of theory [J . This is a testament to the good qual- 
ity of the Burkatski et al. pseudopotentials and ba- 
sis sets employed in our calculations. The geometries we 
used from [f| are sufficiently close to those of Ref. Q that 
they should not appreciably affect the energy differences. 
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TABLE II: The convergence of the QMC energies for structures (a) and (b) with respect to the threshold r\ used to select the 
CSFs included in the trial wave function. The determinant part of the wave function comes from the MCSCF/CAS(8,8) wave 
function with a double-zeta basis. Statistical errors on the last digit are shown in parentheses. 
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TABLE III: DMC energies for different Trotter time steps St, 
using a CAS trial wave function and a double-zeta basis. To- 
tal energies, E a ,E b , of structures (a) and (b) are in Hartrees, 
and the energy difference, AE — E b — E a , is in kcal/mole. 
Statistical errors of the last digit are shown between paren- 
thesis. AE has a negligible time-step error because the errors 
in E a and E b nearly cancel. 
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TABLE IV: Energy differences, AE = E b -E a (in kcal/mole) 
obtained from BLYP, B3LYP, MCSCF/CAS(8,8) methods, 
and, from VMC and DMC using wave functions constructed 
from the BLYP, B3LYP, MCSCF/CAS(8,8) methods. For 
MCSCF both double-zeta (DZ) and triple-zeta (TZ) bases are 
used. The geometries are the same as in Table 1. Statistical 
errors of the last digit are shown in parentheses. 
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Tabic IIVI shows also the QMC energies obtained using 
DFT and CAS (8, 8) determinants to construct the QMC 
wave functions. Note that the QMC energy differences 
obtained using the two single determinants in the trial 
wave functions agree within statistical error, despite the 
fact that B3LYP favors structure (b) and BLYP favors 
structure (a). It appears that it is necessary for this 
system to employ multi-determinant wave functions to 
obtain the correct energy ordering in QMC. 

The best DMC result is obtained with the CAS(8,8) 
wave function; this wave function captures the static cor- 
relations that are important for this system and has a 



lower energy and lower root-mean-square fluctuations of 
the energy than the single determinant wave functions. 
We show in Tablc lTVl thc values obtained with the doublc- 
and triple-zeta basis sets; there is no dependence on the 
basis size within statistical error. Our best estimate is 
that structure (a) has a lower energy than structure (b) 
by 1.9(2) kcal/mole. This result is in excellent agree- 
ment with the best existing quantum chemistry calcu- 
lations, namely 1.0 kcal/mole from CCSD(T) using a 
6-311++G(3df,3pd) basis [| and w 1.6 kcal/mole from 
CCSD(T) using a cc-pVTZ basis @. The result is also in 
reasonable agreement with the 4.7 kcal/mole energy dif- 
ference obtained using a CASSCF(8,8)/6-31+G(2df,2p) 
wave function followed by second order perturbation the- 
ory (CASPT2) 0. 

In summary, our DMC results support the CCSD(T) 
and CCSDT conclusions [J, [f| that the cyclic structure 
(a) is the ground state of the system. We find that accu- 
rate fixed-node DMC energy differences for the structures 
of m-benzyne require multi-determinant wave functions; 
the fixed-node error in the energy difference with a sin- 
gle determinant from restricted B3LYP or BLYP calcu- 
lations is rs 5 kcal/mole. Our best DMC estimate, ob- 
tained with a reoptimized MCSCF/CAS(8,8) trial wave 
function, shows that the proposed monocyclic structure 
of m-benzyne is lower in energy than the bicyclic struc- 
ture by 1-9(2) kcal/mole, in very good agreement with 
previous coupled cluster estimates. 
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